Cours_6

Je prépare ma session de travail

setwd("~/GitHub/UNIGE/32M7129/Cours_06")

Je charge les packages utiles

if(!require("stylo")){
  install.packages("stylo")
  library(stylo)
}

## Loading required package: stylo

## 
## ### stylo version: 0.6.9 ###
## 
## If you plan to cite this software (please do!), use the following reference:
##     Eder, M., Rybicki, J. and Kestemont, M. (2016). Stylometry with R:
##     a package for computational text analysis. R Journal 8(1): 107-121.
##     <https://journal.r-project.org/archive/2016/RJ-2016-007/index.html>
## 
## To get full BibTeX entry, type: citation("stylo")

if(!require("cluster")){
  install.packages("cluster")
  library(cluster)
}

## Loading required package: cluster

if(!require("FactoMineR")){
  install.packages("FactoMineR")
  library(FactoMineR)
}

## Loading required package: FactoMineR

if(!require("factoextra")){
  install.packages("factoextra")
  library(factoextra)
}

## Loading required package: factoextra

## Loading required package: ggplot2

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

if(!require("smacof")){
  install.packages("smacof")
  library(smacof)
}

## Loading required package: smacof

## Loading required package: plotrix

## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

## 
## Attaching package: 'smacof'

## The following object is masked from 'package:base':
## 
##     transform

Ma première CAH depuis 0

Une classification ascendante hiérarchique.

# Je crée un tableau avec des données au hasard
tableau <- cbind(x=c(1,2,5,8,11),y=c(2,4,1,12,10))
colnames(tableau)<- c("Avoir", "Être")
rownames(tableau) <- c("Racine","Molière","Scudéry","Scarron","Corneille")
tableau

##           Avoir Être
## Racine        1    2
## Molière       2    4
## Scudéry       5    1
## Scarron       8   12
## Corneille    11   10

Calcul de distance, par défaut euclidienne. On peut utiliser le paramètre method=“manhattan” pour une autre distance

distance_tableau <- dist(tableau)
distance_tableau

##              Racine   Molière   Scudéry   Scarron
## Molière    2.236068                              
## Scudéry    4.123106  4.242641                    
## Scarron   12.206556 10.000000 11.401754          
## Corneille 12.806248 10.816654 10.816654  3.605551

Classification/clustering. Par défaut complete linkage. On peut utiliser le paramètre method=“centroid” ou “ward.D”

hc <- hclust(distance_tableau)
plot(hc)

#On customise le graph
plot(hc, main="Ma première classification",  xlab="Euclidian distance, \n Complete-linkage clustering" , sub="" , hang =-1)

On peut faire soi-même un delta de Burrows. Pour cela je z-score mes données

#Je transforme mon tableau en data.frame
tableau_df<-as.data.frame(tableau)
#Je z-score mes deux colonnes avec la fonction scale()
tableau_z_a<-scale(tableau_df$Avoir)
tableau_z_b<-scale(tableau_df$Être)
#Je reforme un tableau, avec les nouvelles données
tableau_z<-cbind(tableau_z_a,tableau_z_b)
tableau_z

##             [,1]       [,2]
## [1,] -1.05786348 -0.7724598
## [2,] -0.81743996 -0.3659020
## [3,] -0.09616941 -0.9757388
## [4,]  0.62510115  1.2603292
## [5,]  1.34637170  0.8537714

Et je produis un nouveau graph

hc <- hclust(distance_tableau, method="ward.D")
plot(hc, main="Ma seconde classification",  xlab="Burrows Delta, \n Ward's method" , sub="" , hang =-1)

Un premier cas

Je charge mon corpus

BOYER_AMOURSJUPITERSEMELE_1666<-paste(scan("corpus/BOYER_AMOURSJUPITERSEMELE_1666.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
BOYER_ARISTODEME_1648<-paste(scan("corpus/BOYER_ARISTODEME_1648.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
CORNEILLEP_ANDROMEDE_1651<-paste(scan("corpus/CORNEILLEP_ANDROMEDE_1651.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
CORNEILLEP_MEDEE_1682<-paste(scan("corpus/CORNEILLEP_MEDEE_1682.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
DURYER_DYNAMIS_1653<-paste(scan("corpus/DURYER_DYNAMIS_1653.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
DURYER_CLITOPHON<-paste(scan("corpus/DURYER_CLITOPHON_1632.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
MOLIERE_AMPHITRYON_1668<-paste(scan("corpus/MOLIERE_AMPHITRYON_1668.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
MOLIERE_MISANTHROPE_1667<-paste(scan("corpus/MOLIERE_MISANTHROPE_1667.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
RACINE_PHEDRE_1697<-paste(scan("corpus/RACINE_PHEDRE_1697.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
RACINE_ESTHER_1689<-paste(scan("corpus/RACINE_ESTHER_1689.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
SCARRON_DOMJAPHETDARMENIE_1653<-paste(scan("corpus/SCARRON_DOMJAPHETDARMENIE_1653.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
SCARRON_ECOLIERDESALAMANQUE_1655<-paste(scan("corpus/SCARRON_ECOLIERDESALAMANQUE_1655.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
ROTROU_HERCULEMOURANT_1636<-paste(scan("corpus/ROTROU_HERCULEMOURANT_1636.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
ROTROU_DOMBERNARDDECABRERE_1647<-paste(scan("corpus/ROTROU_DOMBERNARDDECABRERE_1647.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
SCUDERY_MORTDECESAR_1637<-paste(scan("corpus/SCUDERY_MORTDECESAR_1637.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
SCUDERY_ORANTE_1636<-paste(scan("corpus/SCUDERY_ORANTE_1636.txt", what="character", sep="", fileEncoding="UTF-8"),collapse=" ")
#Je peux regarder ici un exemple (enlever le # au début de la ligne)
#CORNEILLEP_ANDROMEDE_1651

J’en fais un tableau

#Je crée une liste de mes textes
my.corpus.raw = list(BOYER_AMOURSJUPITERSEMELE_1666,
                     BOYER_ARISTODEME_1648,
                     DURYER_DYNAMIS_1653,
                     DURYER_CLITOPHON,
                     MOLIERE_AMPHITRYON_1668,
                     MOLIERE_MISANTHROPE_1667,
                     SCARRON_DOMJAPHETDARMENIE_1653,
                     SCARRON_ECOLIERDESALAMANQUE_1655)
#Je tokenise mon texte (espace, apostrophe, retrait de la ponctuation)
my.corpus.clean = lapply(my.corpus.raw, txt.to.words)
#Je compte la fréquence des tokens
complete.word.list = make.frequency.list(my.corpus.clean)
#Je transforme le résultat en table, où les résultats de chaque mot/texte sont alignés
table.of.frequencies=make.table.of.frequencies(my.corpus.clean, complete.word.list, relative = F)

## processing  8  text samples

## 

## combining frequencies into a table...

#Je donne un nom à mes colonnes
row.names(table.of.frequencies)=c("boyer_agamemnon",
                                  "boyer_amoursJupiter",
                                  "duryer_dynamis",
                                  "duryer_clitophon",
                                  "moliere_amphytrion",
                                  "moliere_misanthrope",
                                  "scarron_domJaphet",
                                  "scarron_ecolierSalamanque")
#Je sauvegarde une copie
write.csv(table.of.frequencies, file = "table.of.frequencies.csv",row.names=TRUE)
#je convertis mes données en dataframe
table.of.frequencies = as.data.frame(read.csv(file="table.of.frequencies.csv", sep = ",", header = TRUE, row.names=1, quote = '\"'))
#J'affiche le résultat
View(table.of.frequencies)

Ce résultat ne m’arrange pas: je vais avoir besoin d’avoir les occurrences en rang, et les textes en colonne, et donc d’inverser mon tableau

#j'utilise la fonction transpose(?)
table.of.frequencies<-t(table.of.frequencies)
View(table.of.frequencies)

J’observe la distribution de mon corpus.

summary(table.of.frequencies)

##  boyer_agamemnon   boyer_amoursJupiter duryer_dynamis    duryer_clitophon
##  Min.   :  0.000   Min.   :  0.000     Min.   :  0.000   Min.   :  0.00  
##  1st Qu.:  0.000   1st Qu.:  0.000     1st Qu.:  0.000   1st Qu.:  0.00  
##  Median :  0.000   Median :  0.000     Median :  0.000   Median :  0.00  
##  Mean   :  2.208   Mean   :  1.425     Mean   :  2.127   Mean   :  1.92  
##  3rd Qu.:  1.000   3rd Qu.:  0.000     3rd Qu.:  1.000   3rd Qu.:  1.00  
##  Max.   :599.000   Max.   :321.000     Max.   :589.000   Max.   :537.00  
##  moliere_amphytrion moliere_misanthrope scarron_domJaphet
##  Min.   :  0.000    Min.   :  0.000     Min.   :  0.000  
##  1st Qu.:  0.000    1st Qu.:  0.000     1st Qu.:  0.000  
##  Median :  0.000    Median :  0.000     Median :  0.000  
##  Mean   :  1.783    Mean   :  2.018     Mean   :  1.686  
##  3rd Qu.:  1.000    3rd Qu.:  1.000     3rd Qu.:  1.000  
##  Max.   :478.000    Max.   :567.000     Max.   :425.000  
##  scarron_ecolierSalamanque
##  Min.   :  0.000          
##  1st Qu.:  0.000          
##  Median :  0.000          
##  Mean   :  1.828          
##  3rd Qu.:  1.000          
##  Max.   :423.000

Je peux synthétiser ce résultat avec une boîte à moustache des fréquences (absolues) des mots du corpus.

100% center

Source: stat4decision

maBoite <- boxplot(colSums(table.of.frequencies), main = "Distribution du nombre de mots par de texte", ylab="Nombre de mots (fréq. absolues)", sub="Corpus: distribution")

Je contrôle l’hétérogénéité de mon corpus

#Je dois regrouper les données par auteur
BOYER = colSums(table.of.frequencies[,grepl('boyer_', colnames(table.of.frequencies))])
DURYER = colSums(table.of.frequencies[,grepl('duryer_', colnames(table.of.frequencies))])
MOLIERE = colSums(table.of.frequencies[,grepl('moliere_', colnames(table.of.frequencies))])
SCARRON = colSums(table.of.frequencies[,grepl('scarron_', colnames(table.of.frequencies))])
#Et le graphique
boxplot(list(BOYER,MOLIERE,DURYER,SCARRON), names=c('BOYER','MOLIERE','DURYER','SCARRON'), main="Longueur des textes", sub="Corpus: distribution",ylab="Nombre de mots (fréquences absolues)")

Je transforme cette table des fréquences absolues en une table des fréquences relatives

#Je fais une copie de ma table de fréquences
freqs_rel = table.of.frequencies
#Dans cette copie, pour chacun des mots de chaque colonne, je divise le chiffre trouvé par la somme de la colonne
for(i in 1:ncol(freqs_rel)){
        freqs_rel[,i] = freqs_rel[,i]/sum(freqs_rel[,i])
}
head(freqs_rel)

##      boyer_agamemnon boyer_amoursJupiter duryer_dynamis duryer_clitophon
## de        0.03211280          0.02666556     0.02459929       0.03310319
## et        0.02728784          0.02475494     0.03278050       0.02404143
## je        0.01613681          0.01918923     0.01775378       0.01658242
## vous      0.02299898          0.01844160     0.01703028       0.00961657
## que       0.01704820          0.01910616     0.02276269       0.02416471
## le        0.01715542          0.01702941     0.02075913       0.01713722
##      moliere_amphytrion moliere_misanthrope scarron_domJaphet
## de           0.03173339          0.03162403        0.02984131
## et           0.02701985          0.03150669        0.02260918
## je           0.02243909          0.02434874        0.02569864
## vous         0.01566753          0.03326684        0.02155596
## que          0.02077939          0.02018306        0.01530684
## le           0.01427339          0.01513729        0.01987080
##      scarron_ecolierSalamanque
## de                  0.02738573
## et                  0.02732099
## je                  0.02596141
## vous                0.02537874
## que                 0.02026415
## le                  0.01573223

Je peux (et même je dois) sélectionner dans cette liste les n plus fréquents

#Je prends les premiers rangs. Changez ce chiffre pour changer les résultats par la suite
freqs_rel_mfw = freqs_rel[1:100,]

Analyse

La CAH

Regardons à quoi ressemble notre CAH

CAH = agnes(as.dist(dist.wurzburg(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
plot(CAH)

Nettoyons tout cela…

plot(CAH,
     main="Cluster analysis",
     xlab=paste("Cosine Delta\n Agglomerative coefficient =", CAH_orig_aggloCoeff),
     hang = -0.1)

## Warning in plot.window(xlim, ylim, log = log, ...): "hang" is not a
## graphical parameter

## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "hang" is not a graphical parameter

## Warning in axis(1, at = at.vals, labels = lab.vals, ...): "hang" is not a
## graphical parameter

## Warning in axis(ax$side, at = 0:(length(x$order) - 1), las = 1, labels =
## labels, : "hang" is not a graphical parameter

Faisons un peu mieux

#on va colorier les vecteurs en fonction de l'auteur (en pratique, la chaîne de caractère avant l'_underscore_)
labels_color = vector(length = length(CAH$order.lab))
labels_color[grep("boyer", CAH$order.lab)] = "darkblue"
labels_color[grep("duryer", CAH$order.lab)] = "darkgrey"
labels_color[grep("scarron", CAH$order.lab)] = "darkred"
labels_color[grep("moliere", CAH$order.lab)] = "darkgreen"

#je reprends mon graph
plot(CAH)

#Je fais un peu de coloriage
factoextra::fviz_dend(CAH,
                      ylab=paste("Cosine Delta\n Agglomerative coefficient =", round(CAH_orig_aggloCoeff, digit=3)),
                      k=2, # je cherche à dégager deux clusters
                      rect = TRUE,
                      k_colors = rep("black",2), #je garde la couleur des traits en noir
                      labels_track_height = 0.7, # Je gère la taille de mon rectangle
                      label_cols = labels_color, #j'applique mes couleurs
                      horiz=T) # je mets mon graph à l'horizontal, parce que c'est plus pratique

C’est bien mieux, mais pas encore parfait: le code est est trop long, pas très propre. Cachons tout cela dans un coin. Pour cela, compliquons-nous le travail un peu plus en créant une fonction je que conserve dans un autre fichier appelé functions.R.

Je vais donc charger cette fonction avant de l’utiliser.

source("functions.R")

Quatre paramètres ont été prévus * Le nom des données * La légende à afficher * Le coeffichient d’agglomération (récupéré au passage, un peu plus haut) * une variable pour le paramètre labels_track_height

customPlot(CAH,"Original texts, 100 MFW, Culled @ 0%, Distance: wurzburg",CAH_orig_aggloCoeff,0.7)

Je peux donc reproduire ma CAH simplement avec d’autres configurations. Tentons désormais avec plus de MFW

#on passe de 100 à 1000
freqs_rel_mfw = freqs_rel[1:1000,]
#je recalcule ma CAH
CAH = agnes(as.dist(dist.wurzburg(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
#Et désormais, je n'ai plus qu'une simple fonction "maison" pour me produire ma visualisation
customPlot(CAH,"Original texts, 1000 MFW, Culled @ 0%, Distance: wurzburg",CAH_orig_aggloCoeff,0.5)

Tentons avec une autre distance: Ruzicka measure ou Minmax (cf.Koppel, M. and Winter, Y. (2014). Determining if two documents are written by the same author. Journal of the Association for Information Science and Technology, 65: 178–87) dont la robustesse a été démontrée par Kestmont et al. (Kestemont, M.Stover, J.Koppel, M.Karsdorp, K.Daelemans, W., Authorship Verification with the Ruzicka Metric, DH 2016, Jul 2016, Krakow (Poland)).

100% center

(tf pour term frequency)

#je recalcule ma CAH
CAH = agnes(as.dist(dist.minmax(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
#Et désormais, je n'ai plus qu'une simple fonction "maison" pour me produire ma visualisation
customPlot(CAH,"Original texts, 1000 MFW, Culled @ 0%, Distance: Minmax",CAH_orig_aggloCoeff,-0.2)

Et encore une autre, en diminuant de nouveau les MFW

#on repasse de 1000 à 100
freqs_rel_mfw = freqs_rel[1:100,]
#je recalcule ma CAH
CAH = agnes(as.dist(dist.wurzburg(t(freqs_rel_mfw))), method="ward")
CAH_orig_aggloCoeff <- CAH$ac
customPlot(CAH,"Original texts, 100 MFW, Culled @ 0%, Distance: wurzburg",CAH_orig_aggloCoeff,0.7)

Under the hood: les classes

Tout cela est très bien. Mais qu’est-ce qui est caché derrière ces clusters? Nous pouvons le savoir. Prenons notre dernière CAH, et regardons le processus de construction des clusters:

CAH_2 = as.hclust(CAH)
plot(CAH_2$height, type="h", ylab="height")

Je remarque que j’ai plusieurs clusters. Tentons de voir ce qu’ils contiennent. Classons par 4 notre corpus:

#Je "coupe en deux" ma CAH, ou plutôt je demande une classification en deux des différents
classes = cutree(CAH, k = "4")
#J'ajoute mes classes aux fréquences que j'ai utilisées
Processed_classes = t(freqs_rel_mfw)
Processed_classes = cbind(as.data.frame(Processed_classes), as.factor(classes))
colnames(Processed_classes)[101] = "Classes"

#J'affiche le résultat
Processed_classes[101]

##                           Classes
## boyer_agamemnon                 1
## boyer_amoursJupiter             1
## duryer_dynamis                  2
## duryer_clitophon                2
## moliere_amphytrion              3
## moliere_misanthrope             3
## scarron_domJaphet               4
## scarron_ecolierSalamanque       4

Diminuons à deux le nombre de cluster

#Je "coupe en deux" ma CAH, ou plutôt je demande une classification en deux des différents
classes = cutree(CAH, k = "2")
#J'ajoute mes classes aux fréquences que j'ai utilisées
Processed_classes = t(freqs_rel_mfw)
Processed_classes = cbind(as.data.frame(Processed_classes), as.factor(classes))
colnames(Processed_classes)[101] = "Classes"

#J'affiche le résultat
Processed_classes[101]

##                           Classes
## boyer_agamemnon                 1
## boyer_amoursJupiter             1
## duryer_dynamis                  1
## duryer_clitophon                1
## moliere_amphytrion              2
## moliere_misanthrope             2
## scarron_domJaphet               2
## scarron_ecolierSalamanque       2

Je peux aussi regarder les valeurs associées à chacun des mots, et ceux associés à chaque cluster:

myClasses = catdes(Processed_classes, num.var = 101)
myClasses

## 
## Link between the cluster variable and the quantitative variables
## ================================================================
##            Eta2      P-value
## je    0.8855821 0.0004896757
## amour 0.8040975 0.0025455982
## yeux  0.7601916 0.0047643342
## moi   0.7512318 0.0053410502
## y     0.7445267 0.0058031317
## pas   0.7375033 0.0063162127
## être  0.7371581 0.0063422162
## n     0.7356576 0.0064561179
## si    0.7193183 0.0077904004
## l     0.7138624 0.0082757444
## ne    0.6959504 0.0100182228
## faire 0.6770794 0.0121175393
## dieux 0.6744536 0.0124323540
## qui   0.6573420 0.0146275157
## mes   0.6342385 0.0180111247
## en    0.6330403 0.0182005359
## veux  0.6181419 0.0206765720
## m     0.6095886 0.0222026475
## ses   0.5960834 0.0247753466
## plus  0.5890029 0.0262070694
## suis  0.5741345 0.0294079634
## comme 0.5457179 0.0363077270
## bien  0.5423877 0.0371876525
## me    0.5098090 0.0466527660
## point 0.5027205 0.0489296469
## 
## Description of each cluster by quantitative variables
## =====================================================
## $`1`
##          v.test Mean in category Overall mean sd in category   Overall sd
## amour  2.372484      0.004818126  0.003380804   0.0008299353 0.0016028752
## yeux   2.306803      0.001808802  0.001422762   0.0001256128 0.0004427618
## si     2.243931      0.009360696  0.007472486   0.0015781275 0.0022263300
## l      2.235405      0.017151679  0.015865132   0.0009675451 0.0015227141
## dieux  2.172826      0.003528194  0.001978623   0.0014225437 0.0018868420
## qui    2.145086      0.009527600  0.008502894   0.0005062741 0.0012638733
## mes    2.107052      0.003316553  0.002528891   0.0007548847 0.0009890388
## ses    2.042690      0.002753466  0.001970091   0.0007580942 0.0010146507
## plus   2.030522      0.005931663  0.005243258   0.0006308718 0.0008969846
## suis  -2.004730      0.001103213  0.002013334   0.0005036370 0.0012011368
## m     -2.065701      0.004543850  0.005795718   0.0008642927 0.0016033934
## veux  -2.080143      0.001166075  0.001530247   0.0001085254 0.0004631936
## en    -2.105061      0.010588157  0.012344905   0.0015522742 0.0022079724
## faire -2.177052      0.001783577  0.002343016   0.0004299235 0.0006798804
## ne    -2.207182      0.008047953  0.009833099   0.0010665234 0.0021398561
## n     -2.269274      0.005291597  0.006619240   0.0008650303 0.0015479011
## être  -2.271587      0.001702773  0.002472563   0.0003633443 0.0008965862
## pas   -2.272119      0.004991906  0.006090712   0.0006925025 0.0012794955
## y     -2.282912      0.001141360  0.001708866   0.0003815097 0.0006577039
## moi   -2.293169      0.003940317  0.005481745   0.0003737930 0.0017784282
## je    -2.489794      0.017415563  0.021013767   0.0011821199 0.0038235906
##          p.value
## amour 0.01766891
## yeux  0.02106579
## si    0.02483681
## l     0.02539075
## dieux 0.02979338
## qui   0.03194597
## mes   0.03511304
## ses   0.04108310
## plus  0.04230349
## suis  0.04499194
## m     0.03885672
## veux  0.03751245
## en    0.03528599
## faire 0.02947669
## ne    0.02730134
## n     0.02325169
## être  0.02311147
## pas   0.02307934
## y     0.02243556
## moi   0.02183829
## je    0.01278171
## 
## $`2`
##          v.test Mean in category Overall mean sd in category   Overall sd
## je     2.489794     0.0246119713  0.021013767   0.0013957565 0.0038235906
## moi    2.293169     0.0070231732  0.005481745   0.0011974515 0.0017784282
## y      2.282912     0.0022763722  0.001708866   0.0002747232 0.0006577039
## pas    2.272119     0.0071895170  0.006090712   0.0006163696 0.0012794955
## être   2.271587     0.0032423534  0.002472563   0.0005390367 0.0008965862
## n      2.269274     0.0079468818  0.006619240   0.0007200349 0.0015479011
## ne     2.207182     0.0116182445  0.009833099   0.0012833566 0.0021398561
## faire  2.177052     0.0029024541  0.002343016   0.0003371910 0.0006798804
## en     2.105061     0.0141016520  0.012344905   0.0010809285 0.0022079724
## veux   2.080143     0.0018944192  0.001530247   0.0003899695 0.0004631936
## m      2.065701     0.0070475865  0.005795718   0.0011226719 0.0016033934
## suis   2.004730     0.0029234554  0.002013334   0.0009875056 0.0012011368
## plus  -2.030522     0.0045548539  0.005243258   0.0005131879 0.0008969846
## ses   -2.042690     0.0011867150  0.001970091   0.0005069224 0.0010146507
## mes   -2.107052     0.0017412296  0.002528891   0.0003817371 0.0009890388
## qui   -2.145086     0.0074781878  0.008502894   0.0009156383 0.0012638733
## dieux -2.172826     0.0004290517  0.001978623   0.0005425611 0.0018868420
## l     -2.235405     0.0145785853  0.015865132   0.0006251135 0.0015227141
## si    -2.243931     0.0055842772  0.007472486   0.0005403116 0.0022263300
## yeux  -2.306803     0.0010367229  0.001422762   0.0002797223 0.0004427618
## amour -2.372484     0.0019434819  0.003380804   0.0005637698 0.0016028752
##          p.value
## je    0.01278171
## moi   0.02183829
## y     0.02243556
## pas   0.02307934
## être  0.02311147
## n     0.02325169
## ne    0.02730134
## faire 0.02947669
## en    0.03528599
## veux  0.03751245
## m     0.03885672
## suis  0.04499194
## plus  0.04230349
## ses   0.04108310
## mes   0.03511304
## qui   0.03194597
## dieux 0.02979338
## l     0.02539075
## si    0.02483681
## yeux  0.02106579
## amour 0.01766891

Under the hood: l’es classes’ACP

Je commence par faire une ACP.

ACP = PCA(t(freqs_rel_mfw))

Faisons un peu de ménage, ajoutons de l’information, et tentons de comprendre

#categories
get_categories = vector(length = length(colnames(freqs_rel_mfw)))
get_categories[grep("racine", colnames(freqs_rel_mfw))] = "Racine"
get_categories[grep("moliere", colnames(freqs_rel_mfw))] = "Moliere"
get_categories[grep("duryer", colnames(freqs_rel_mfw))] = "Duryer"
get_categories[grep("boyer", colnames(freqs_rel_mfw))] = "Boyer"

ACP = PCA(t(freqs_rel_mfw))

fviz_pca_ind(ACP, col.ind = get_categories, legend="none")

On voit nettement deux axes, qu doivent bien s’appuyer sur quelque chose. Serait-il possible de faire ressortir les mots qui sont à l’origine de ce partitionnement des données?

Continuons notre nettoyage pour y voir plus clair

#categories
get_categories = vector(length = length(colnames(freqs_rel_mfw)))
get_categories[grep("racine", colnames(freqs_rel_mfw))] = "Racine"
get_categories[grep("moliere", colnames(freqs_rel_mfw))] = "Moliere"

ACP = PCA(t(freqs_rel_mfw))

fviz_pca_ind(ACP, col.ind = get_categories, legend="none")

fviz_pca_var(ACP, col.var="contrib",geom.var = "text", select.var = list(contrib =20))+
scale_color_gradient2(
                      low="yellow", mid="orange", 
                      high="darkred", midpoint=1.6)+theme_bw()

Nous nous rappelons qu’il existe un pourcentage de réalité retenu par axe, que nous pouvons observer, et donc calculer. Dans ce cas précis, il y en a seulement trois axes, car nous n’avons mis que quatre de textes.

ACP$eig

##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1  30.037372              30.037372                          30.03737
## comp 2  21.041028              21.041028                          51.07840
## comp 3  15.055366              15.055366                          66.13377
## comp 4  11.101088              11.101088                          77.23485
## comp 5   9.201767               9.201767                          86.43662
## comp 6   7.811868               7.811868                          94.24849
## comp 7   5.751511               5.751511                         100.00000

Un petit graph montre bien la lente perte de significativité des axes :

barplot(ACP$eig[,1], main="Percentage of variance", names.arg=1:nrow(ACP$eig))

Les autres algorithmes

t-SNE

Avec l’algorithme de Barnes-Hut,

library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

Attention! Si on relance la même commande, on obtient un résultat différent!

library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#encore une fois
library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#encore une fois
library(Rtsne)
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 0.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

On peut faire varier la perpléxité (c’est à dire s’autoriser une distorsion plus grande): on passe à 1.5

maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 1.5, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

On peut faire varier la perpléxité: on passe à 2

maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

#Encore une fois
maRtsne = Rtsne(t(freqs_rel_mfw), dims = 2, initial_dims = 36, perplexity = 2, theta = 0.0, check_duplicates = TRUE, pca = TRUE)
plot(maRtsne$Y)
text(maRtsne$Y[,1], maRtsne$Y[,2], labels = row.names(t(freqs_rel_mfw)), cex=.6) 

MDS

Le Multidimensional scaling est une compression des données en deux dimensions. Il en existe trois types: * classique * métrique * non métrique

MDS classique

Le MDS classique se base sur le calcul de distance, dont il va essayer de préserver l’essentiel.

#Je crée mes données
fit = cmdscale(dist(t(freqs_rel_mfw), method = "manhattan"), eig=TRUE, k=2) # k est le nombre de dimensions souhaité
#Je dessine mon résultat
x = fit$points[,1]
y = fit$points[,2]
plot(x, y, xlab=paste("Coordinate 1 (GOF: ", round(fit$GOF[1] * 100, digits=2), "%)"), ylab=paste("Coordinate 2 (GOF: ", round(fit$GOF[2] * 100, digits=2), "%)"), main="PMD métrique")
text(x, y, labels = row.names(t(freqs_rel_mfw)), cex=.7) 

MDS Métrique

Le MDS métrique, dit aussi ordinal, ne s’intéresse pas à la mesure de distance, mais sa valeur en relation avec les autres paires

Rappelons que le stress permet de mesurer la distortion introduite dans le résultat.

#J'affiche le stress avec le paramètre verbose=T(RUE)
MDSmetrique = mds(dist(t(freqs_rel_mfw), method = "manhattan"), ndim=2, type="ordinal",verbose=T)

## Iteration:    1  Stress (raw):   0.01510422  Difference:    0.06233185 
## Iteration:    2  Stress (raw):   0.00887573  Difference:    0.00622848 
## Iteration:    3  Stress (raw):   0.00656184  Difference:    0.00231389 
## Iteration:    4  Stress (raw):   0.00541200  Difference:    0.00114984 
## Iteration:    5  Stress (raw):   0.00475869  Difference:    0.00065331 
## Iteration:    6  Stress (raw):   0.00439233  Difference:    0.00036636 
## Iteration:    7  Stress (raw):   0.00417426  Difference:    0.00021807 
## Iteration:    8  Stress (raw):   0.00403217  Difference:    0.00014209 
## Iteration:    9  Stress (raw):   0.00393218  Difference:    0.00009999 
## Iteration:   10  Stress (raw):   0.00385734  Difference:    0.00007485 
## Iteration:   11  Stress (raw):   0.00379958  Difference:    0.00005776 
## Iteration:   12  Stress (raw):   0.00375696  Difference:    0.00004262 
## Iteration:   13  Stress (raw):   0.00372424  Difference:    0.00003271 
## Iteration:   14  Stress (raw):   0.00369790  Difference:    0.00002634 
## Iteration:   15  Stress (raw):   0.00367617  Difference:    0.00002173 
## Iteration:   16  Stress (raw):   0.00365798  Difference:    0.00001819 
## Iteration:   17  Stress (raw):   0.00364262  Difference:    0.00001537 
## Iteration:   18  Stress (raw):   0.00362954  Difference:    0.00001308 
## Iteration:   19  Stress (raw):   0.00361836  Difference:    0.00001119 
## Iteration:   20  Stress (raw):   0.00360875  Difference:    0.00000961 
## Iteration:   21  Stress (raw):   0.00360047  Difference:    0.00000828 
## Iteration:   22  Stress (raw):   0.00359331  Difference:    0.00000716 
## Iteration:   23  Stress (raw):   0.00358711  Difference:    0.00000620 
## Iteration:   24  Stress (raw):   0.00358173  Difference:    0.00000538 
## Iteration:   25  Stress (raw):   0.00357706  Difference:    0.00000467 
## Iteration:   26  Stress (raw):   0.00357299  Difference:    0.00000407 
## Iteration:   27  Stress (raw):   0.00356944  Difference:    0.00000355 
## Iteration:   28  Stress (raw):   0.00356633  Difference:    0.00000311 
## Iteration:   29  Stress (raw):   0.00356361  Difference:    0.00000272 
## Iteration:   30  Stress (raw):   0.00356122  Difference:    0.00000238 
## Iteration:   31  Stress (raw):   0.00355913  Difference:    0.00000209 
## Iteration:   32  Stress (raw):   0.00355730  Difference:    0.00000184 
## Iteration:   33  Stress (raw):   0.00355568  Difference:    0.00000161 
## Iteration:   34  Stress (raw):   0.00355427  Difference:    0.00000142 
## Iteration:   35  Stress (raw):   0.00355302  Difference:    0.00000125 
## Iteration:   36  Stress (raw):   0.00355193  Difference:    0.00000110 
## Iteration:   37  Stress (raw):   0.00355096  Difference:    0.00000096

plot(MDSmetrique, sub=paste("Stress, ", round(MDSmetrique$stress, digits=2)))

Non métrique

#J'affiche le stress avec le paramètre verbose=T(RUE)
MDSnonMetrique = mds(dist(t(freqs_rel_mfw), method = "euclid"), ndim=2, type="interval",verbose=T)

## Iteration:    1  Stress (raw):   0.02349141  Difference:    0.05037561 
## Iteration:    2  Stress (raw):   0.01787966  Difference:    0.00561176 
## Iteration:    3  Stress (raw):   0.01617650  Difference:    0.00170315 
## Iteration:    4  Stress (raw):   0.01541740  Difference:    0.00075910 
## Iteration:    5  Stress (raw):   0.01499855  Difference:    0.00041885 
## Iteration:    6  Stress (raw):   0.01474414  Difference:    0.00025441 
## Iteration:    7  Stress (raw):   0.01458330  Difference:    0.00016085 
## Iteration:    8  Stress (raw):   0.01447937  Difference:    0.00010393 
## Iteration:    9  Stress (raw):   0.01441109  Difference:    0.00006828 
## Iteration:   10  Stress (raw):   0.01436553  Difference:    0.00004556 
## Iteration:   11  Stress (raw):   0.01433468  Difference:    0.00003085 
## Iteration:   12  Stress (raw):   0.01431351  Difference:    0.00002117 
## Iteration:   13  Stress (raw):   0.01429880  Difference:    0.00001471 
## Iteration:   14  Stress (raw):   0.01428848  Difference:    0.00001033 
## Iteration:   15  Stress (raw):   0.01428115  Difference:    0.00000732 
## Iteration:   16  Stress (raw):   0.01427593  Difference:    0.00000523 
## Iteration:   17  Stress (raw):   0.01427217  Difference:    0.00000376 
## Iteration:   18  Stress (raw):   0.01426945  Difference:    0.00000272 
## Iteration:   19  Stress (raw):   0.01426747  Difference:    0.00000197 
## Iteration:   20  Stress (raw):   0.01426603  Difference:    0.00000144 
## Iteration:   21  Stress (raw):   0.01426498  Difference:    0.00000105 
## Iteration:   22  Stress (raw):   0.01426421  Difference:    0.00000077

plot(MDSnonMetrique, sub=paste("Stress, ", round(MDSmetrique$stress, digits=2)))